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Abstract 

A parallelized hybrid Monte Carlo (HMC) methodology is devised to quantify the microstructural 
evolution of poly crystalline material under elastic loading. The approach combines a time explicit 
material point method (MPM) for the mechanical stresses with a calibrated Monte Carlo (cMC) 
model for grain boundary kinetics. The computed elastic stress generates an additional driving 
force for grain boundary migration. The paradigm is developed, tested, and subsequently used to 
quantify the effect of elastic stress on the evolution of texture in nickel polycrystals. As expected, 
elastic loading favors grains which appear softer with respect to the loading direction. The rate 
of texture evolution is also quantified, and an internal variable rate equation is constructed which 
predicts the time evolution of the distribution of orientations. 
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I. INTRODUCTION 



The properties of polycrystalline materials, ranging from chemical resistance to fracture 
toughness, depend on the orientations and populations of the grains of which they are com- 
prised jl L 2]. When these orientations are non-random, the material is said to have material 
texture 3(. This anisotropic feature of materials often develops as a result of manufacturing 
and associated heat treatment [lj]. Texture can also result from cyclic thermomechanical 
loading which causes some orientations to grow at the expense of others. Such cyclic load- 



ing is common, for instance, in most electrica 



on the solder connections of circuit boards 



devices where it can have a deleterious effect 
5[. Thermally induced loading leads to elastic 



stresses which can be partially relieved by the growth of preferentially aligned grains. This 
type of microstructural evolution, at relatively high temperatures and within the realm of 
very small deformations, is referred to as stress-induced texture. This texture has not been 
studied to the same degree as the texture development associated with materials process- 
ing 3], 6-8], and it is the focus of the work presented here. Localized yield may also occur, 
of course, but such plastic deformation is not considered in the current study. 

Meso-scale computational interrogations of texture evolution are now commonplace {9- 
We are particularly interested, though, in grain boundary motion in which elastic driving 
forces are not negligible. Even within this narrow focus, phase-field 9| and sharp-interface 
(SI) (l2| approaches have been explored. Missing from the arsenal, though, are Monte Carlo 
(MC) methods, and this is most likely because a discrete, probabilistic setting does not lend 
itself to the calculation of stress. The current work remedies this by developing a hybrid 
algorithm in which a deterministic, continuum computation of elastic energy is dovetailed 
with MC kinetics for microstructural evolution. In order to focus on the methodology, 
we restrict our attention to evolution dominated by grain boundary kinetics and disregard 
independent triple junction kinetics {13I. 1^ ]. This implies that the grain boundary properties 
are isotropic. We further assume that the boundary mobility is not a function of its accretive 
velocity, and so the von Neumann conditions are met [3]. Simulations carried out with a 



high MC temperature, which is different 
naturally satisfy these requirements [ill . 



rom the temperature of the physical material, most 
3- 



The physical domain is idealized as a continuum with sharp, coherent interfaces with 
material properties quantized on a cubic grid. MC algorithms generate probabilities for grid 



point changes in crystal orientation based upon the associated change in global free energy. 
A quasi-static elastic contribution to the free energy change can therefore be generated by 
comparing stress states before and after each trial event. We present a numerically efficient 
means of approximating elastic contribution by merely calculating the stress in a small zone 
around the trial position. The boundaries of the zone are given constant traction boundary 
conditions and the work done by the zone on the rest of the domain is accounted for in the 



16]. 



free energy shift. This is implemented for a Potts model 

A previously developed calibration procedure, here extended from two to three dimen- 
sions, is used to endow the hybrid Monte Carlo (HMC) paradigm with physical length, 



time, and energy scales based on experimentally measurable SI properties [15(]. In order 
to efficiently implement the approach within a parallelized environment, the computational 
domain is decomposed into a checkerboard and orientation is updated using a Red/Black 
(RB) algorithm [17j. An alternative would be to adopt the N-Fold Way 18|, |l9| updat- 
ing strategy which would dramatically speed up a serial MC calculation but has yet to be 
parallelized j^ . 

The HMC method is applied to investigate the microstructural evolution of nickel poly- 
crystals under elastic loading. Although our primary task is to develop the meso-scale 
paradigm, we also show that the data generated can be used to construct a macroscopic 
kinetic equation for the evolution of the variance in a Gaussian distribution of grain orien- 
tations. This equation can be used to efficiently estimate the effects of elastic loading on 
material texture. 

II. GRAIN BOUNDARY KINETICS WITHIN THE ELASTIC REGIME 

SI theory is briefly reviewed because it is used to endow the MC paradigm with physical 
time and length scales. We focus on a polycrystalline cube of linear elastic, anisotropic 
material. Away from grain boundaries, the elastic distortion energy, U, stress S, and strain 
e are related: 

U=^S-e, S = Ce. (1) 
The fourth order elastic stiffness tensor, C, depends on grain orientation. Under our idealized 



assumptions, the thermodynamic driving traction for interfacial accretion is 



12|: 



f s = [[£/]] - <S> ■ [[e]] + k 7 * . (2) 



Here 7* is the capillary driving force, k is the local mean curvature. The expressions <~> 
and [["]] represent domain average and difference across the grain boundary, respectively. 
The Herring relation 2l| can then be used to describe the accretive normal speed of the 
interface, v: 

v = m s f s . (3) 

The proportionality constant, m s , is grain boundary mobility, which can be measured ex- 
perimentally along with driving force The information is then translated into a discrete, 
probabilistic setting for the prediction of grain boundary motion. Specifically, Equation ([3]) 
is replaced by a fluctuation based flip probability which is carried out at all discrete points 
on grain borders. 



III. HYBRID MONTE CARLO METHODOLOGY 



In the absence of elastic stresses, texture development is assumed to be captured by a 



Q-state Potts model 16| introduced here within a non-dimensional context. The unit cube 
domain is divided into a uniform (K x K x K) lattice. The non-dimensional reference lattice 
spacing is A = 1/K. Grain orientation is described with an integer- valued spin field, 
which ranges from 1 to Q. One MC time step is defined as the completion of N trial events 
(N = K 3 ). There is no time scale inherent in the MC model, so it is assigned a time interval 
of r mc for use in subsequent calibration. The system Hamiltonian is assumed to be 



N M N 
i=l n=l i=l 

where 6j is the spatially varying elastic energy associated with the i th lattice, J quqn is the 
interaction energy between the two neighbor bins, M is the number of neighbors considered 
of the selected lattice, and 5 is the Kronecker delta. 

Kinetics, within the MC paradigm, is treated via a series of trial changes in the orientation 
of randomly selected sites. The probability with which such trial flips are accepted is 
related to the overall change in the Hamiltonian of the system. It will be convenient to 
use the following notation for the change in a field before and after a flip in the phase (spin 
variable) of one site: [[H]]* = H Trml — H. Here H Trml is the Hamiltonian after the flip event. 



A standard Metropolis algorithm 



22| is adopted wherein the probability, P, for each flip is 



a function of the resulting change in Hamiltonian: 




[[%]]* >0 

M* < o 



(5) 



The fundamental MC temperature, T mc , parameterizes variability and is not the physical 
temperature of the material. An issue is how to account for the elastic contribution to [[%]]*■ 
Suppose that a flip probability is to be calculated at a single discrete point on a grain 
boundary within a domain of infinite extent. The elastic contribution to [[H]]* is simply 
the change in the total elastic energy, U, integrated over the entire domain. In principle, 
this must be calculated by equilibrating the stress twice — once before the flip and once 
afterwards. This is computationally intractable, and so an approximation is made. First 
create a neighborhood, fix, around the point and consider the rest of the domain to be a 
reservoir, f2 2 (Fig. d]). A trial flip is assumed to not affect the traction on the boundary 
of the interior region. Under this assumption, the equilibrated elastic energy which results 
from a trial flip is the sum of the elastic energy of Q± and the work done by Hi to expand 
or contract its boundary under constant traction. This boundary motion results from a 
trial change of elastic stiffness associated with the new orientation being considered. This 
relation approaches the exact trial energy when Qi increases to include the entire domain. 
The elastic contribution to [[H]]* is therefore approximated as 



where Tj and Ui represent traction and displacement associated with the trial flip, respec- 
tively. 

The most extreme case is to reduce the sub-domain Q± to a single point so that 



where p and i represent the selected particle and its trial state, respectively. Then we have 




(6) 




(7) 




(8) 



IV. TIME, LENGTH AND ENERGY SCALES 



In order to consider the elastic driving force in MC simulations, a correspondence must be 
established between the parameters of the MC and SI paradigms. This correspondence has 

5 



Conservative 
Reservoir 

fi 2 



Material 



Applied Forces — ► 
(need not be same for each site) 



FIG. 1. The material domain and its conservative reservoir. 



been developed for a generalized bulk energy in 2D 15j, and here we extend the approach 
to a 3D, elastic setting. The characteristic energy, Eq, length, Iq, and time, to are needed to 
derive expressions for experimentally measurable grain boundary stiffness, 7%, bulk energy 
density, b s , and grain boundary mobility, m s , in terms of their Monte Carlo counterparts, 
7*, b m and m: 



1*E 
A 2 J/2 
- _ 2b m E 

The MC time unit, r mc , and interaction energy, J, are typically set to unity. It is then 
assumed that both experimental and MC parameters are specified, and the expressions are 
used to solve the characteristic energy, length, and time scales which link the SI paradigm 
to a physically meaningful SI setting. A key field in these links is the MC grain boundary 
mobility, m(a,tp), which is calculated using 

m{aM = ^ L y (10) 

The parameter, m*(a), is the MC grain boundary reduced mobility, and it can be measured 
by studying the capillarity driven shrinkage of internal grains 

of energy, the nondimensional parameter, a = ^p-, characterizes variability within the 
MC system. The physical domain size is taken to be the characteristic length, l Q . The 
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characteristic time, to is then, 

fc=S=£^, (11) 

and the characteristic energy, E , is written as 

^* I 2 A 2 

Eo = ^4^- (12) 

Here 7% is the experimentally measured grain boundary stiffness, and the MC grain bound- 
ary stiffness, 7*, can be numerically derived 151 ] . 

The SI elastic driving force, [[K]]*, is then converted into a MC setting using 

mw m = B^5. (13) 

Here is the elastic driving force in MC format. This expression for elastic driving 

force is then used in Eq. flSJ) to give a probabilistic rate of grain boundary motion strictly in 
terms of SI parameters. This endows the MC paradigm with length, time and energy scales. 



V. ELASTIC DEFORMATION 



A time-exp 



mat ion 



23 



icit version of Material Point Method (MPM) is used to track elastic defor 



251 ] . Material is viewed 



24| , although other algorithms could be complementary 
as a set of discrete mass points, and a back ground mesh is set to solve the momentum 
equation. In response to a prescribed loading condition, the mechanical simulator breaks 
boundary loading into a series of sub-loadings. For each sub-loading, the time-explicit algo- 
rithm is launched to seek the quasi-static state in the following manner. Nodal velocity is 
computed and mapped back to material points. The strain increment associated with each 
material point is obtained by calculating the gradient of its velocity. Hooke's law is then 
applied to update stresses of each material point. The state variables of material points are 
subsequently re-mapped to nodes to start the next iteration. This procedure continues until 
either the relative stress increment or strain increment criterion is met. The relative stress 
increment, Act, and the strain increment, Ae, of each particle are defined as 

ct* - a'- At 

and 

^ = e\- e\-^\ (15) 
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FIG. 2. The dependence of the Ni Young's modulus on the rotation angle with respect to the 



z-axis. 



respectively. Here % represents the particle index, t is the current time step and (t — At) 
is the previous time step. The quasi-static state is reached when either the relative stress 
increment or relative strain increment meets the prescribed convergence criterion for all the 
particles at time t. 

Attention is restricted to cubic symmetry for which there are only three independent 
elastic constants, Cn, Cyi and C44. Grain orientation at each computational grid point 
determines the principal frame for the elastic stiffness, C: 

C = Cijkiei <g> e,j <g> e fc <g> e h (16) 

Here {e^} is an arbitrary Cartesian basis, and the element Cyfej can be written by the index 



relationship between matrix and tensor expressions 26]. A rotated stiffness tensor, C, is 
then given by 

^■'ijkl = / ^y^y^y^ ^--tifc^imAj n Afc A;p, (17) 

m n o p 

where A is a general rotation matrix, which is written as 

A = BCD. (18) 
The individual rotation matrices, B, C, and D can be parameterized by three Euler angles, 

As an example, subsequently considered in more detail, a nickel single crystal is deformed 
under uniaxial loading with the elastic moduli of Cn = 250.8, Cu = 150.0 and C44 = 123.5 
GPa. The effective Young's modulus along the x-axis of the rotated crystal is depicted in 
Fig. [2] as a function of crystal rotation about its z-axis. 
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The static poly crystalline Ni response to elastic loading was then computed. To ensure 



isotropic behavior, the number of grains, N„, was determined using 



2n 



N 9 =°-^( V -1), (19) 



^iso 



where e; so is the acceptable relative error of calculated stress and strain among independent 
simulations, and rj is the characteristic parameter of anisotropy: 

2a_jAA , . 

V = r r ■ (20) 

^11 — <^12 

This implies that isotropy along with a 2% error will be achieved for a nickel polycrystal 
with 300 grains. The experimental value of the Young's modulus of polycrystalline Ni is 225 
GPa {29)], while our MPM calculation gives 221 GPa. The difference is 1.8%. The agreement 
between simulations and experiments allows us to apply the technique to a polycrystalline 
setting in which the grain boundaries evolve. 

VI. MICROSTRUCTURAL EVOLUTION 

In typical implementations of the Q-state Potts simulations, updating follows immediately 



after a trial event is accepted. This amounts to a Gauss-Seidel (GS) update [la]. The 
efficiency of this approach can be remarkably improved using the N-Fold Way algorithm 19 ]. 
When applied within a parallel environment, though, communication between processors will 
be asynchronous. Since we seek to develop a simulation approach amenable to large-scale, 
parallelized environments, we give up the N-Fold Way in favor of a Red/Black (RB) updating 
scheme fl7 |. In the RB scheme, an MC domain is decomposed into a checkerboard. Since 
each spin site interacts only with its four nearest neighbors, the outcome of a trial event 
is determined only by the four opposite colored lattice sites. The update of the selected 
particle will not affect any other particle with a different color. This allows the particles 
updating in parallel. One MC step is therefore split into black/red and red/black half steps. 
In each, the trial states of the red or black particles are not updated until this half step 
finishes. Only two updating steps are required in RB algorithm as opposed to a series of 
updates associated with GS algorithm. Orientations on the boundaries of the decomposed 
domain require a communication between processors. 

The RB updating a ppr oach is validated by recovering the predictions of Isotropic Grain 
Growth Theory [j], 0,130, 3l| which predicts that the mean grain radius at time t is R = 
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FIG. 3. Relationships between rate of change of average grain size and MC time for two updating 
algorithms. The results pertain to a 3D polycrystal. The solid gray line and the solid black lines 
are fits to the GS and RB data, respectively. Grid size=60 3 and a = 1.0. Error bars are for the 
standard deviation of the 10 averaged simulation results. Volume is measured as the number of 
cubic grid cells comprising a grain. 




cit 1/n . Here c\ is a constant and n is referred to as the grain growth exponent. Within the 
current setting, n — 2 Therefore, grain volume, V, is predicted to evolve as V 3 ^ 2 = c 2 t 1 
where c 2 is a constant. 

A set of three-dimensional simulations were considered on a 60 3 grid, wherein a total of 
99 orientations were randomly distributed. Grain boundaries were all given the same inter- 
action energy, and the microstructure was subsequently allowed to evolve with a variability 
parameter a = 1.0. The average grain size was collected over time and is shown in Figj3l 
The data indicate that both RB and GS have the correct power law behavior, although the 
mobility of the RB algorithm is higher than that for GS. This is because the first half step 
updating in RB algorithm creates more boundary corrugation and the capillarity effect is 
therefore accelerated. 



VII. RESULTS 

The HMC algorithm was first analyzed within an idealized one-dimensional setting. A 
bi- crystalline bar of length, L, and cross-sectional area, A, was fixed at left and loaded with 
constant traction, F, at right. The Young's moduli of the two domains at left and right are 
Ei and E 2 (Ei > E 2 ), respectively. The domain is discretized by 100 lattice sites (N = 100). 
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Eq. (|S]) was then used to calculate the elastic driving force: 

M =± £<i-i>- (21) 

Positive driving force is obtained if a particle of grain two is selected and flipped to grain 
one, and the trial flip is therefore accepted. This implies the grain boundary moves to the 
left side with a speed of one particle distance per MC step, i.e., 

L 



(22) 



Nt 

The HMC implementation delivered a grain boundary velocity (averaged over 100 simula- 
tions) which is consistent with this analytical solution, and the relative error was found to 
be less than one percent. Note that the grain boundary moves so as to grow the softer 
grain. This is the guiding rule for texture evolution in response to elastic loading; softer 
orientations will grow at the expense of their stiffer counterparts. Here stiffer means that 
the grain offers a higher effective Young's modulus along the axis of loading. Of course, 
both elastic stored energy and work are considered. This mechanism is different from the 



approaches which consider the elastic stored energy only [32]. 

A series of three-dimensional simulations were subsequently carried out. One mm 3 poly- 
crystalline nickel cube was deformed at 450 °C. For computational convenience, in later 
calibration, grain boundary energy and grain boundary mobility are taken to be isotropic, 
although in reality the two items are characterized by anisotropy [l|, Q]- Room temperature 
elastic moduli were adopted even though there is a 8% difference between room temperature 



and 450 °C values |29|]. One million computational particles were used with each assigned 
a random initial orientation. Only the rotation angle with respect to the z-direction was 
varied to more clearly interpret the results. The other two Euler angles were set to zero. The 
rotation angle was varied from 0° to 90° in 1° steps; therefore 91 orientations were used in 
total. Fig. [2] shows the orientation effect on the Young's modulus in the x-direction. The 45° 
is the stiffest orientation, while 0° and 90° are the softest orientations. The random initial 
microstructure was first allowed to coarsen over 100 MC steps at a = 2.6 to construct a 
polycrystalline microstructure. A uniaxial tension of 30 MPa was applied in the x-direction 
in 75 loading substeps. In each sub-step, the quasi-static elastic state was obtained and was 
then followed by one MC step during which the microstructure evolved. In general, tex- 
ture evolution is described by two mechanisms, grain boundary kinetics or triple junction 
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Parameters 


Values 


* 

m 


0.27 


7* 


2.4 


m s (m 4 /J.s) 


1.0 x 10" 13 


7% (J/m 2 ) 


0.45 


m* (m 2 /s) 


0.45 x 10~ 13 


l (m) 


1.0 x 10~ 3 


to (a) 


600 


£ (J) 


1.875 x 10" 11 



TABLE I. Grain boundary properties and required characteristic parameters to link MC and SI 
models. 

kinetics |14|. In order to more cleanly isolate the influence of elastic loading, we focus only 
on the former and thus implicitly assume that grain boundary properties are isotropic. As 
previously noted, we also assume that the the mobility is not a function of accretive veloc- 



ity M . Within the MC setting, these conditions are met by restricting attention to high 
MC temperatures. The constant, experimentally reasonable, grain boundary energy [1|, 133[ 



and grain boundary mobility 34J are assigned, although the energy and mobility are known 
to depend on the details of grain boundary structure [l[ . Key boundary properties and their 
corresponding MC parameters, obtained from Eq. ( fTTi) and Eq. ( fT2l) . are listed in Table HI 

Since any sufficiently high MC temperature is valid, we chose a value such that one MC 
step is equal to 10 minutes. Fig. |4] shows the Young's modulus probability density associated 
with individual orientations at three time slices. Consistent with intuition and the results 
of our ID study, evolution favors grains which are softer with respect to the elastic loading 
axis. To quantify the elastic effect on texture evolution, the texture histograms were fitted 
to a time dependent Gaussian equation: 



i 



g(t)V^Erf(45/(V2g(t))) 



Exp[- 



(23) 



where / is the orientation probability density, <p is l ne orientation angle, and Erf is an 
error function. The evolution of texture is thus distilled to a single variance, g. In order 
to characterize this process with a single parameter, the rate of change of the variance is 
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FIG. 4. Evolution of texture under elastic loading with a single Euler angle description of each 
orientation. Every five of 91 simulation data are shown as discrete points along with Gaussian fits 
as solid curves. Results are averaged over 10 independent runs. Errors bars quantify the standard 
deviation. 



assumed to be proportional to the current variance, 

where r is a characteristic time constant. Since Equation (12 ip indicates that the driving 
force is proportional to a 2 , the solution to Equation (124"]) is written as 

g{t) = Po*e-$, (25) 

where the time constant, r, and coefficient, /3, were then fitted to the HMC data: r = 
500 min and p = 0.05 °/MPa 2 . Within our restrictive ansatz, the characteristic time, r, 
describes the rate of texture evolution in polycrystalline Ni under a uni-axial loading. This 
fitted Gaussian variance, g(t), was subsequently used in Eq. fT23|) to produce the solid curves 
of Fig. 1 

More general polycrystalline settings can not use a single angle of misorientation to 
describe texture, so we now show a second approach which employs the effective Young's 
modulus as a texture measure. First the Young's modulus data of Fig. [2] are fitted to a 
Gaussian function: 

(6 - 45°) 2 

x{<t>) = aExp[- W J ] + c, (26) 

where x((f>) = E(<f>)/C 11 , a = 0.388, b = 16.40°, and c = 0.541. It is the inverse of this 
relation, though, that is required in order to describe the texture in terms of the Young's 
modulus: 

0( x ) = 45° - {-2b 2 Log(— -)]*, 0.55 < x < 0.929. (27) 

(X 

13 



10 

200 300 400 500 600 700 
Time (Min) 

FIG. 5. Evolution of the Gaussian variance for the distribution of orientations. The solid dots are 
the HMC results while the gray curve is the solution of Eq. (|24p . 
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FIG. 6. Re-plot of texture data of Fig. U]with horizontal axis now the effective Young's modulus. 



The orientation probability density, f, can now be written as a function of Young's modulus 
by combining Equations (127]) and (1231) . and the result is plotted in Fig. [61 The distribution, 
now referred to as the Young's modulus probability density, is initially uniform but favors 
softer grains with increasing time. 

The previous numerical experiment used only one Euler angle in order to make the study 
particularly transparent. A more general application was therefore considered in which 999 
random orientations span all three Euler angles. Fig. [7] shows the resulting texture evolution 
quantified by the effective Young's modulus. The data are fitted to the composition of 
functions given in Equations (12 7p and (1231) with a = 0.667, b = 16.40°, and c = 0.541. Note 
that only the parameter a differs from its value for a single angle description of orientation, 
and the previously obtained Gaussian variance, g, is utilized. 
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FIG. 7. Evolution of texture under elastic loading with three arbitrary Euler angles description of 
each orientation. The 999 data points are integrated into 91 points, and every three of the 91 data 
points are shown as discrete points along with Gaussian fits as solid curves. Results are averaged 
over 20 independent runs. Errors bars quantify the standard deviation. 

VIII. DISCUSSION 

The hybrid Monte Carlo (HMC) paradigm combines a deterministic, continuum elasticity 
solver with a probabilistic, discrete algorithm for grain boundary evolution. The paradigm 
is intended for application where elastic loading influences grain boundary motion on time 
scales sufficiently slow that inertial effects can be neglected. 

The isothermal, uniaxial loading of polycrystalline Ni is adopted in to demonstrate the 
idea. Attention was restricted to a MC temperature regime for which the grain boundary 
properties are isotropic and variations in grain orientation were initially restricted to single 
Euler angle. In line with intuition, it was found that grains which are softer with respect 
to the loading direction are energetically favored. The results also show that the texture 
conforms to a Gaussian histogram with a variation that can be described by a first order 
differential equation. For quantifying the texture evolution where more general Euler angles 
are allowed, the independent variable can be changed from angle to the effective Young's 
modulus with respect to loading direction. 

The HMC methodology was then applied to a more realistic Ni polycrystal in which all 
three Euler angles were used. As expected, the softer grains with respect to the loading 
direction eventually dominate the polycrystalline media. The data were accurately fitted 
to the same relation used for a single angle of orientation. The fit describes the rate at 
which texture evolves using a single characteristic time. Significantly, characteristic time 
at which the distribution tightens was found to be the same as for the case where only a 
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single angle of orientation was allowed to vary. An appropriate next step is to compare 
the predictions of this model directly with experimental data for texture evolution under 
uniaxial loading. We are not aware of any such data with which we can currently make 
such a comparison. In all implementations, only the most primitive approximation was used 
for the elastic contribution to energy changes do to re-orientation. The domain over which 
stress equilibration is considered could be extended to a larger neighborhood of trial flip 
sites. 

This hybrid computational methodology for studying texture evolution offers an alter- 
native to sharp-interface and phase-field modeling. The resulting texture evolution maps 
are expected to be useful in materials and process design as well as in part performance 
assessment throughout service life. An extension of this work to include linearized inelastic 
deformation, currently underway, will allow the influence of linear plastic deformation on 
texture to be accounted for as well. The approach can also be extended to account for an 
evolving temperature field via coupling with the heat equation. 
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X. APPENDIX. NUMERICAL ERROR ANALYSIS 



The numerical errors associated with MC simulations have been comprehensively stud- 



ied 22], so here we focus on errors associated with the mechanical response 



mbalances 



in force can be made arbitrarily small by tightening the convergence criteria 35]. In prac- 
tice, though, limits should be set which balance the desire to reduce such numerical errors 
with the need to obtain an efficient algorithm. The two interior convergence parameters 
of relative stress increment and strain increment were set to 1.0 x 10~ 9 and 1.0 x 10 -5 , 
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FIG. 8. Uniaxial tension test results of a 2D polycrystal: (a) computed engineering strains and 
(b) the relative errors of the engineering strains with respect to grid size. 



respectively. A second class of numerical errors is associated with the background mesh size 
and the number of computational particles per cell. To quantify these, uniaxial tension tests 
were carried out on a polycrystalline system with seven columnar, hexagonal grains. The 
total number of computational particles was varied but the number of particles per cell was 
fixed at four, i.e., two particles per cell in each direction in a given cross-section. Figs. M 
(a) and (b) demonstrate the calculated engineering strains and relative errors, respectively. 
The minimum number of particles required to meet a 1% relative error was determined to 
be 12 x 10. This analysis provides a criterion to determine the number of particles required 
to represent one grain in polycrystalline media. A direct extension to 3D implies that 900 
computational particles should be used per grain, and this rule is followed in all simulations. 



XI. REFERENCES 



[1] F. J. Humphreys and M. Hatherly. Recrystallization and Related Annealing Phenomena. 
Elsevier Ltd, Oxford, United Kingdom, 2004. 

[2] U. F. Kocks, C. N. Tome, and H.-R. Wenk. Texture and Anisotropy. Cambridge University 
Press, Cambridge, United Kingdom, 1998. 

[3] I. L. Dillamore and W. T. Roberts. Rolling textures in f.c.c. and b.c.c. metals. Acta Metal- 
lurgical 12(3) :281 - 293, 1964. 



[4] A. U. Telang, T. R. Bieler, A. Zamiri, and F. Pourboghrat. Incremental recrystallization/grain 
growth driven by elastic strain energy release in a thermomechanically fatigued lead- free solder 
joint. Acta Materialia, 55(7):2265 - 2277, 2007. 

[5] G. Hofer. Texture Dependent Young's Modulus in Austenitic Cladding. Textures and Mi- 
cro structures, 8-9:611-617, 1987. 

[6] J. Hirsch and K. Lucke. Overview no. 76: Mechanism of deformation and development of 
rolling textures in polycrystalline f.c.c. metals-i. description of rolling texture development in 
homogeneous cuzn alloys. Acta Metallurgica, 36(11):2863 - 2882, 1988. 

[7] K. Lucke, R. Rixen, and M. Senna. Formation of recrystallization textures in rolled aluminum 
single crystals. Acta Metallurgica, 24(2):103 - 110, 1976. 

[8] G. B. Sarma, B. Radhakrishnan, and T. Zacharia. Finite element simulations of cold defor- 
mation at the mesoscale. Computational Materials Science, 12(2):105 - 123, 1998. 

[9] L. Q. Chen. Phase-Field Models for Microstructure Evolution. Annual Review of Material 
Research, 32:113-140, 2002. 
[10] M. E. Gurtin and M. T. Lusk. Sharp-interface and phase-field theories of recrystallization in 

the plane. Physica D: Nonlinear Phenomena, 130(1-2):133 - 154, 1999. 
[11] E. A. Holm, G. N. Hassold, and M. A. Miodownik. On misorientation distribution evolution 

during anisotropic grain growth. Acta Materialia, 49(15):2981 - 2991, 2001. 
[12] R. Abeyaratne and J. K. Knowles. On the driving traction acting on a surface of strain 
discontinuity in a continuum. Journal of the Mechanics and Physics of Solids, 38(3) :345 - 
360, 1990. 

[13] G. Gottstein, Y. Ma, and L.S. Shvindlerman. Triple junction motion and grain microstructure 

evolution. Acta Materialia, 53(5): 1535 - 1544, 2005. 
[14] G. Gottstein and L.S. Shvindlerman. Grain boundary junction engineering. Scripta Materialia, 

54(6):1065 - 1070, 2006. 

[15] P. Liu and M. T. Lusk. Parametric links among Monte Carlo, phase-field, and sharp-interface 

models of interfacial motion. Phys. Rev. E, 66:061603, 2002. 
[16] F. Y. Wu. The Potts model. Reviews of Modern Physics, 54(l):235-268, 1982. 
[17] H. Fried. The checkerboard update Glauber model, cellular automata and Ising models. J. 

Phys. A, 23:4165-4181, 1990. 



18 



[18] A. B. Bortz, M. H. Kalos, and J. L. Lebowitz. A new algorithm for monte carlo simulation of 

ising spin systems. Journal of Computational Physics, 17(1): 10 - 18, 1975. 
[19] G. N. Hassold and E. A. Holm. A fast serial algorithm for the finite temperature quenched 

Potts model. Computers in Physics, 7(1):97-107, 1993. 
[20] G. Korniss, M. A. Novotny, and P. A. Rikvold. Parallelization of a dynamic monte carlo 

algorithm: A partially rejection- free conservative approach. Journal of Computational Physics, 

153(2):488 - 508, 1999. 

[21] C. Herring. Surface tension as a motivation for sintering. McGraw-Hill, New York, 1949. 

[22] D. Landau and K. Binder. A Guide to Monte Carlo Simulations in Statistical Physics. Cam- 
bridge University Press, Cambridge, United Kingdom, 2005. 

[23] D. Sulsky and H. L. Schreyer. Axisymmetric form of the material point method with applica- 
tions to upsetting and taylor impact problems. Computer Methods in Applied Mechanics and 
Engineering, 139(l-4):409 - 429, 1996. 

[24] A. R. York II. Development of Modifications to The Material Point Method for the Simulation 
of Thin Membranes, Compressible Fluids and Their Interactions. PhD thesis, Sandia National 
Laboratories, 1997. 

[25] A. Needleman, R. J. Asaro, J. Lemonds, and D. Peirce. Finite Element Analysis of Crystalline 
Solids. Computer Methods in Applied Mechanics and Engineering, 52:689-708, 1985. 

[26] T. C. T. Ting. Anisotropic Elasticity: Theory and Applications. Oxford University Press, 
United States, 1996. 

[27] G. Arfken. Mathematical Methods for Physicists. Academic Press, Orlando, FL, 1985. 

[28] M. Nygards. Number of grains necessary to homogenize elastic materials with cubic symmetry. 
Mechanics of Materials, 35:1049-1057, 2003. 

[29] H. M. Ledbetter. Temperature behaviour of Young's moduli of forty engineering alloys. Cryo- 
genics, pages 653-656, 1982. 

[30] E. A. Holm, D. J. Srolovitz, and J. W. Cahn. Microstructual evolution in two-dimensional 
two-phase polycrystals. Acta Materialia, 41(4):1119-1136, 1993. 

[31] V. Tikare, E. A. Holm, D. Fan, and L.Q. Chen. Comparison of phase-field and Potts models 
for coarsening processes. Acta Materialia, 47(1):363-371, 1998. 

[32] B. Radhakrishnan and G. Sarm. Coupled Finite Element-Potts Model Simulations of Grain 
Growth in Copper Interconnects. Mater. Res. Soc. Symp. Proc, 1156-1161, 2009. 



[33] D. L. Olmsted, S. M. Foiles, and E. A. Holm. Survey of computed grain boundary properties 
in face-centered cubic metals: I. grain boundary energy. Acta Materialia, 57(13) :3694 - 3703, 
2009. 

[34] R. Le Gall, G. Liao, and G. Saindrenan. Experimental Determinaton of Nickel Grain Boundary 
Mobility During Recrystallization. Materials Science Forum, 294-296:509-512, 1999. 

[35] J. B. Scarborough. Numerical Mathematical Analysis. Oxford University Express, Oxford, 
London, 1955. 



20 



